Spin Versus Charge Density Wave Order in Graphene-like Systems 
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A variational technique is used to study sublattice symmetry breaking by strong on-site and 
nearest neighbor interactions in graphene. When interactions are strong enough to break sublattice 
symmetry, and with relative strengths characteristic of graphene, a charge density wave Mott in- 
sulator is favored over the spin density wave condensates. In the spin density wave condensate we 
find that introduction of a staggered on-site energy (quasiparticle mass) leads to a splitting of the 
fermi velocities and mass gaps of the quasiparticle spin states. 
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The possibility of gapping the spectrum of graphene, 
either by expHcit^^ or spontaneous sublattice symmetry 
breaking^., is an important fundamental and practical 
problem'^. At the fundamental level, the question of 
spontaneous breaking of either exact or approximate chi- 
ral symmetry emulates similar issues in quantum field 
theories such as quantum chromodynamics. At the prac- 
tical level, a small gap, particularly one which could 
be switched on and off would be important for using 
graphene in electronics technology as it could give a 
mechanism for controlling the flow of electrons. 

In the absence of magnetic fields, the best clean, sus- 
pended graphene is a semi-metal with no discernible 
energy gap. Gap formation by spontaneous symme- 
try breaking, if it occurred, would be driven by strong 
electron-electron interactions. Numerical Monte Carlo 
computations and series expansions of the Hubbard 
model on a hexagonal lattice indicate that a phase transi- 
tion from a semi-metal to an anti-ferromagnetic, or spin 
density wave (SDW), Mott insulator— (with perhaps 
other exotic phases in between) will occur for relatively 
strong coupling, U/t ~3-5, where U is the on-site Hub- 
bard interaction and t is the hopping parameter. Esti- 
mates of these parameters for suspended graphene, where 
an on-site Coulomb energy is [/ ~ lOeV and t = 2.7eV 
have a ratio in the same range, raising the tantalizing 
idea that graphene is close to this critical point and 
some small modification which enhances the interaction 
could induce a phase transition to a gapped siaXe^^—. 
The simplest gapped states are spin density wave and 
charge density wave (CDW) Mott insulators, although 
more exotic phases have been discussed^iii"— . There is 
also a possibility of breaking the sublattice symmetry ex- 
plicitly by depositing graphene on the appropriate sub- 
strate, such as boron nitride or silicon carbidelj^ii^. Even 
once it is broken explicitly, there can be phase transition 
between different patterns, for example CDW to SDW, 
which can be of great interest. Moreover, the interplay 
between spontaneous and explicit symmetry breaking is 
an interesting problem which has been discussed in recent 
literaturei^^i^. 

In this Letter, we shall show that, even with explicit 
symmetry breaking, electron-electron interactions can 




FIG. 1: The phase diagram of the extended Hubbard model 
with staggered potential m. The thick lines are phase bound- 
aries between the SDW phase (sgnA^ = — sgnA^) and the 
SM/CDW phase (sgnA^ = sgnA^), while the thin line for 
m = is the boundary between the SM and the CDW phases. 
The SM phase does not appear when m is finite. As m be- 
comes larger, the SDW phase is suppressed. 



change the character of the gap and the electron spec- 
trum significantly. For example, a candidate for the 
gapped phase is an antiferromagnetic SDW Mott insula- 
tor, and it is indeed what is found in the Hubbard model 
at strong coupling^"— . We shall show that, when next-to- 
nearest neighbor (NN) interactions are added to the Hub- 
bard model, with the strength appropriate to graphene 
{V ^ lOeV), a CDW state is favored over the antifer- 
romagnetic SDW. A quantum phase transition between 
the two can be driven by varying either the strength of 
the NN coupling or the amplitude of an explicit symme- 
try breaking staggered potential. Our central result is the 
phase diagram in Fig.[T] The critical Hubbard coupling is 
underestimated by our technique, likely because magnon 
fluctuations are not taken into account. In the absence of 
explicit symmetry breaking, the semi-metal phase occurs 
in the trapezoid in the lower left-hand corner. 

We shall use a variational technique where we replace 
the full Hamiltonian if by a solvable trial Hamiltonian 
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FIG. 2: The hexagonal graphene lattice composed of sublat- 
tices A (black dots) and B (white dots) connected by the basis 
vectors s^. 



Hq which is optimized using Jensen's inequalitj***. 
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We shall adjust Hq to minimize this upper bound on the 
free energy. Here, (O)o = TTe-^"'>0/TTe-^"" . 

For the Hamiltonian of graphene, we begin with the 
tight-binding model with nearest-neighbor (NN) hop- 
ping. 

HT = -t H{r)ba{r + s,) + bi{r + s,)a^{r)] . 

i,iT,r£A 

(2) 

The hexagonal graphene lattice is depicted in Fig. [21 
It contains two triangular sublattices, A and B. Cre- 
ation and annihilation operators for electrons at sites r on 
sublattice A are (aj^(r), ao-(r)) and B are (6j^(r), &cr(i"))- 
a =ti i is the spin index. We shall add a staggered on-site 
energy, Hm, which models explicit sublattice symmetry 
breaking (which could arise by interaction with a sub- 
strate, for example, and gives the low energy graphene 
Dirac electron a mass gap—), a Hubbard interaction Hjj 
and a NN interaction Hy, 

HM^mJ2 [bU^ + si)ba{r + si) ~ ai{r)aa{r)] (3) 
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(4) 

Hv^vY, [4(i-)a^(r)-l] \bl,{v + s,)b„,{v + s,) - I 



(5) 



where terms such as aJ.(r)ao-(r) are summed over spins. 
An important symmetry of graphene which is to a good 
approximation visible in angle-resolved photoemission 
spectroscopy (ARPES) measurementa^i is particle-hole 
symmetry. Here, we have written a model Hamilto- 
nian H — Ht + Hm + Hjj + Hv which has exact 
particle-hole symmetry. We will also restrict the vari- 
ational Ansatz to have this symmetry. The explicit 



particle-hole transformation is a^(r), ao-(r), &^(r), 6o-(r) 
— ?> ao-(r), aj^(r), —b„{Y),—b\{T). The terms in the Hamil- 
tonian Ht, Hjj, Hv, Hm are invariant. 

To write down the trial Hamiltonian Hq, it is conve- 
nient to Fourier transform to momentum space where 



;(k) -A,(k) ; V ; ■ 

(6) 

where k is a wave- vector in the Brillouin zone of the 
triangular lattice, and, for example 

a(k) = 2^ —^a^{r) , a^{r) = J dk-^a^(k) (7) 

with fl the volume of the Brillouin zone. Here we as- 
sume that the different matrix elements in the Hamilto- 
nian can be simultaneously diagonalized in spin. This 
is not the most general possible Ansatz, which would 
have a more complicated spin dependence. We have as- 
sumed translation invariance on the triangular sublat- 
tices. If we set A„(k) = and h„{k) = ^ e^k s. = 
^(k), Hq becomes identical to the tight-binding model 
Hamiltonian Ht- We have fixed the diagonal parts 
of Hq so that it has particle-hole symmetry. Aside 
from particle- hole symmetry, Ht, Hj/, Hy also have 
sublattice symmetry ~ where we simply interchange 
the sublattice excitations a^(k), acr(k), &J.(k), 6cr(k) — 
6j.(— k), k), aj^(— k), acr(— k). This symmetry is bro- 
ken by Hm, which flips sign under the transformation. 
The trial Hamiltonian has this symmetry only when 
Ao- — 0. Hermiticity requires that /i*(— k) = ^.^(k) and 
A<,(k) = A,(-k) =real. 

The spectrum and the eigenstates of Hq are easy to 
find: The eigenvalues of the single-particle Hamilto- 
nian are E^,±{k) = ±S,(k) - i^Mk^HMkyP- 
With a change of variables into polar coordinate h ~ 
EcoaOe"!' , A = Esin0 {-n/2 < 6* < n/2, - n < (f) < 
tt), Hq is diagonalized by the canonical transformation 
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v/2(l + sin( 
1 



[(1 + sin e)ij+ - cos 6le"^ V-] (8) 
[cosee-"t'ij+ + (1 + sin6l)V'-] ,(9) 



where we have suppressed k, cr labels, and (■0+,'0+) and 
("0- J 4'- ) a-re creation and annihilation operators for elec- 
trons in energy states +i?o-(k) and — i?o-(k), respectively. 
With this transformation, the correlation functions are 
diagonal in momentum and spin space, 

(ai(k)a,(k))o = i [l-sin0,(k)tanh|^,(k)] (10) 

(fet(k)6,(k))o = i [l + sin0,(k)tanh|^,(k)l (11) 



{bl{k)a^{k))Q = -icos6l<,(k)e*'^-(''Hanh 



.(12) 



All the others can be obtained from these by simple al- 
gebra. All expectation values of operators factor into 
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bilinears such as these. Then, the free energy per unit 
volume is the sum of the following five contributions, 
which come from Fq — {Hq)o and the expectation val- 
ues of Ht,Hm,Hu,Hv, respectively: 



tanh^^-lln 



2 cosh 



(13) 

2 J il ^ 

(14) 
(15) 



U 
4 



— E sin 0a (k) tanh f £'<T(k) 
E J ^ sin 0,(k) tanh f£;,(k) 
jE sin f?,(k) tanh |£;,(k: 

— sin^?,(k) tanhfS,(k) 



(16) 
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dk 



cos e„ (k)e''^' $(k) tanh (k) 



(17) 



First, consider the equation obtained from varying (k); 

= Z^e*'^"('')$(k) - Z*e-*"*-(*')$*(k), (18) 
where the factor is defined by 

Za = l-^J ^ cos 0,(k)e-'^^-('')<i>(k) tanh ^i?,(k). 

(19) 

The solution of this equation which minimizes the en- 
ergy is (paik) = — arg$(k) -f tt. Thus, everywhere in 
Eqs. (HSll-dlll), 6'"^$ can be replaced by -|$|. The equa- 
tion obtained by varying £'a(k) and 0a- (k) are 

£;a(k) =cos0,(k)Z„t|$(k)| 

+ [^a-m+i^Ca]sin0a(k), (20) 
Z,i|$(k)| tan0,(k) = - m + ^^C^, (21) 

where 

/dk B 
— sine, (k) tanh |i;,(k) (22) 

^- = 1 + 7^ / ^ cos 0,(k')|<i>(k)| tanh j£;,(k). (23) 
The solution reads 



i?.(k) = V^^t^l$(k)P+[f a-: 

(24) 

sin0,(k) = - m + ^^Cff] /^a(k) (25) 

cos0a(k) = Z,t|$(k)|/i;,(k), (26) 



m/t m/t 

FIG. 3: The behavior of the density wave amplitude Ao- (left) 
and the velocity renormalization factor Za (right) as functions 
of the external mass m, where the on-site interaction U = 6.0t 
and the NN interaction V = 0.5t are fixed. The system shows 
the SDW phase for m < 1.8t, while it reveals the CDW phase 
for m > 1.8t. The fermi velocity of the up spin and that of 
the down spin differ (i.e. ^ Z^) in the SDW phase, unless 
m = 0. 



The four constants Ca and Z^ must be determined self- 
consistently. Zn corrects the fermi velocity and Ca and 
m gap the spectrum. If m were zero, but Ca nonzero, 
the sublattice symmetry would be spontaneously broken. 
The nonzero temperature is important for deriving the 
variational equations, however, to study the low temper- 
ature limit, we will set it to zero. 

Since, from Eq. (ITB|) . eu = ■j^t^'j,, the Hubbard inter- 
action favors a spin density wave (SDW) where C^ and 
Ci are nonzero and have opposite signs. From Eq. pT]) . 
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{C^ + C^Y-2,V[{Z^-lY + {Z^-lY]. TheNN 
interaction favors a charge density wave (CDW) where 
C-f- and Cj, are nonzero and have the same sign. The 
competition of these two phases is seen in the numerical 
solutions of the self-consistent equations, Eqs. ([20|) - ([26| . 
The phase diagram is shown in Fig. [TJ When m = 0, 
there are three phases: a semi-metal (SM) for C/, F < i, 
SDW for U >V,m, and CDW for V,m>U. Hm is a 
source for CDW. When it is finite, there is no SM phase. 
The SDW phase is suppressed as m increases, while the 
CDW phase is enhanced. 

Now we shall investigate the quantitative behavior 
of Act and Za, by varying one parameter out of U, 
V, m while holding the others fixed. We begin with 
U — 6.0t, V = O.bt, m = 0, where the system is in 
the SDW phase, and we increase m. As shown in the left 
panel of Fig. [31 both and A^ increase as a function of 
m as long as m is sufficiently small. It should be noted 
that |A-|-| and \Ai\ take different values in this region 
unless TO = 0, which means that the quasiparticle gap 
for up spin and that for down spin are different. Such 
a discrepancy of |Act| also causes the discrepancy of the 
factor Za through Eq. ([^5]) , as shown in the right panel of 
FigEl The difference between Z^ and Z^ increases as a 



function of to towards its maximum value Z 



i' 



0.05 



at m = 1.7t, then drastically drops towards zero at the 
critical value TOc = I .St. Since A^ = Aj^ in the CDW 
region, quasiparticles with up spin and those with down 
spin obtain the same Fermi velocity above mc- 

Next we vary the NN interaction V, where the on-site 
interaction U = 6.0t and the mass m = O.bt are fixed. 
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FIG. 4: The behavior of the density wave amplitude (left) 
and the velocity renormalization factor Zc (right) as functions 
of the NN interaction V, where the on-site interaction \J = 
6.0t and the external mass m = 0.5t are fixed. The system 
shows the SDW phase for V < 1.5t, while it reveals the CDW 
phase for V > 1.5t. The fermi velocity of the up spin and 
that of the down spin differ (i.e. Zf ^ Z^) in the SDW phase, 
unless V = 0. 

The SDW amplitude is suppressed as V increases, and 
a phase transition to the CDW phase occurs at Vc = 
l.bt, as shown in the left panel of FigU) Due to the 
finite external mass, m, there is a discrepancy between 
|A-|-| and jAj,] in the SDW phase, which leads to the 
discrepancy between and as shown in the right 
panel of FigU) Since Z^ — Ws proportional to V , Z^ and 
Z^ take the identical value (unity) at ^ = 0, even though 
the quasiparticle gap amplitudes are different. — Z^ 
reaches its maximum value 0.09 just below the critical 
value Ve- 
in conclusion, we note that, in the continuum limit of 
graphene, the CDW and SDW condensates are indistin- 
guishable as they are related to each other by a transfor- 
mation in the emergent U(4) symmetry. We have found 
that they are indeed distinguished by lattice scale physics 
which can have an important effect. We have shown 
that the short ranged interactions of relative strengths 
approximating graphene favor the CDW state. This is 
basically due to the fact that the on-site energy is anoma- 
lously small compared to the NN potential energy. Ex- 
plicit symmetry breaking, which can be present in some 



cases enhances this effect. The conclusion that lattice 
scale physics can drive a phase transition is surprising. 
It is likely that the naive continuum Coulomb interac- 
tion is good for the semi-metal phase, however when den- 
sity wave order sets in, it is driven by otherwise irrele- 
vant four-Fermion interactions which can have nontriv- 
ial strong coupling fixed points. This point of view is 
supported by renormalization group analyses of the con- 
tinuum theory^. Another anomalous effect of explicit 
symmetry breaking is the splitting of the Fermi veloci- 
ties of the spin up and spin down electrons in the SDW 
phase. That splitting goes to zero if m goes to zero. It 
increases as a function of the NN interaction strength, 
and it reaches about 10% of the Fermi velocity. Such a 
discrepancy might be detected by ARPES measurements, 
and it may influence transport properties of the system. 

Some of our results are similar to a self-consistent mean 
field theory. We point out that the variational technique 
is more general in that it contains a wave-function renor- 
malization, which is normally absent in mean field ap- 
proach. It is also readily applicable to a much wider 
array of potentials, and we believe that our exposition 
of the technique here could be used as a starting point 
for more general analyses of graphene-like systems. We 
have focused on the SDW and CDW patterns, but the 
honeycomb lattice can have a richer array of symmetry 
breaking patterns, such as the Kekule distortion which is 
expected to become relevant when the next-to-NN inter- 
action is taken into account. The interplay of ordering 
patterns including those phases, induced either sponta- 
neously or explicitly, remains an open question. 
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